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Modern distributed systems include a class of applications in which non-functional requirements 
are important. In particular, these applications include multimedia facilities where real time 
constraints are crucial to their correct functioning. In order to specify such systems it is necessary 
to describe that events occur at times given by probability distributions and stochastic automata 
have emerged as a useful technique by which such systems can be specified and verified. 

However, stochastic descriptions are very general, in particular they allow the use of general 
probability distribution functions, and therefore their verification can be complex. In the last few 
years, model checking has emerged as a useful verification tool for large systems. In this paper 
we describe two model checking algorithms for stochastic automata. These algorithms consider 
how properties written in a simple probabilistic real-time logic can be checked against a given 
stochastic automaton. 
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1. INTRODUCTION 

In this paper we describe and compare two model checking algorithms for stochastic 
automata. The reason for building such model checking algorithms is to support 
the verification of non-functional properties in distributed multimedia systems. 

The advent of distributed multimedia applications such as video conferencing, 
collaborative virtual environments, video on demand etc, place great demands on 
the specification and design of such systems because of the need to describe and 
verify non- functional requirements Bowman ct al. 1996[ ]. These non-functional 



requirements typically involve real time constraints such as placing bounds on end- 
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to-end latency, and are often called Quality of Service (QoS) [Bowman et al. 1998] 
requirements because they reflect the overall quality of delivery as opposed to the 
functional aspects. 

In order to specify and verify such constraints it is necessary not only to be 
able to describe deterministic timing concerns but also probabilistic and stochastic 
systems. That is, in practice timings cannot be assumed to be fixed (determinis- 
tic timings) but events can occur at different times with particular probabilities. 
Therefore it is necessary to describe timings that occur according to certain prob- 
ability distributions. For example, in a network specification it is not sufficient to 
assume that the packet deliveries arrive at fixed predetermined times, instead we 
need to model the system where they might arrive at times determined by (for 
example) an exponential distribution. 

There are now a number of techniques which can be used to describe such sys- 
tems, e.g. Queueing Systems [ Kleinrock 1975||, Ge neralised Stochastic Petri-nets 
[Marsam et al. 1984], Markov Chains [Stewart 1994 ], generalised semi-Markov pro- 
cesses flGlynn 198C Sto chastic Process Algebra jHillston 1996[ and Stochastic Au- 
tomata flD'Argcnio 1999 ] etc. In this paper we consid er Stochastic Automata (which 
are related to timed automata [Alur and Dill 1994 ) . We define two model check- 
ing | Baier and Kwiatkowska 1998 algorithms for them. 

Stochastic automata are a very promising specification and verification paradigm. 
In particular they allow the study of both functional and non-functional require- 
ments within the same description, giving a more complete view of overall perfor- 
mance than, say, a queueing theory description of the problem. They also support 
not just exponential distributions but general distributions. The issue here is the 
following. In a stochastic specification we need to associate a distribution function 
F with an action a so that we can describe the probability of the time delay after 
which a can happen. Stochastic automata naturally allow general distributions, 
in contrast say to stochastic process algebras which usually restrict themselves to 
exponential distributions jHillston 1996]. 

In practice it is unrealistic to only consider exponential distributions and it is 
necessary for arbitrary distributions (e.g. uniform, gamma, deterministic etc) to 
be considered. For example, it is often assumed that packet lengths are exponen- 
tially distributed. However, in reality this is not the case, rather they are cither 
of constant length (as in ATM cells Tanenbaum 1996 1) or they are uniformly dis- 
tributed with minimum and maximum size (as in Ethernet frames [Tanenbaum 



1996]). Stochastic automata allow such arbitrary distributions to be used. 



There are ostensibly two ways to move from the tractable case of exponential 
distributions to the less tractable case of generalised distributions. One approach 
is to make small generalisations of markov chains by allowing limited forms of 
non-memoryless behaviour (see e.g. GSPNs | Marsam et al. 1984 |). However, the 
problem with this approach is that there will always be classes of distributions 
that cannot be modelled. The alternative is to allow any distribution, but then 
use heuristics and coarse approximation techniques to contain the problem of in- 
tractability. The majority of work on this topic follows the first of these approaches. 
However here we investigate the feasibility of the second approach and thus we im- 
pose few constraints on the generality of the distributions we allow in our stochastic 
automata. 

Because stochastic automata are related to timed automata, verification strate- 
gies for stochastic automata can be derived by using the extensive work on verifica- 
tion for timed automata, see e.g. [Larsen et al. 1997] |Daws et al. 1995 1 [Henzinger 
et al. 1997 1 . The p articular verification technique we consider is model check- 
ing |Alur et al. 199Cf ]. This is perhaps the most successful technique to have arisen 



Stochastic Model-Checking for Multimedia 



3 



System 
(Automaton) 



Property 
(Temporal Logic 
Formula) 



Model Checker 



~y 

System models formula 



System does not 
model formula 



Fie. 1. Model checker 



from concurrency theory. The basic approach is to show that an automaton de- 
scription of a system satisfies a temporal logic property, see Figure 



In accordance with a number of other workers, e.g. | Baier et al. 1999fl , we view the 
application of model checking to analysis of stochastic systems as a very exciting 
combination, since it provides a form of generalised transient analysis — for exam- 
ple the property [->error Li<woo error] < 0.01 states that the probability of first 
reaching an error state within 1000 time units is less than 1 percent, and whether 
a particular stochastic system satisfies this property can be investigated. 

In defining our model checking algorithm we draw heavily on the experience of 



model checking timed automata e.g. |Larsen et al. 1997 1. However, the move from 



timed to stochastic leads to new issues that must be tackled. In particular, many 
of the properties that we wish to verify are inherently probabilistic. Conventional 
model checking allows us to answer questions such as "Is a particular sequence of 
events possible?", but in stochastic model checking we want to ask "What is the 
probability of this sequence of events?" . To do this we will check an automaton 
against a simple probabilistic temporal logic. 

We present two approaches to model checking stochastic automata. Both ap- 
proaches are enumerative in the sense that, in showing whether a property holds, 
they enumerate reachable configurations of the system. However, the methods by 
which they determine the probability of being in a particular configuration are quite 
different. Specifically, one derives probabilities by integrating the relevant proba- 
bility density functions, while the second responds to the difficulties incurred in 
evaluating these integrals (which will become clear during the paper) by employing 
a discretisation process. 

The structure of the paper is as follows. In Section |] we introduce stochastic au- 
tomata illustrated by a simple example. In Section]^ we define a small probabilistic 
real-time logic, in which we can express simple properties that we wish to check 
our stochastic automata against. The first algorithm is presented in Section ^, and 
the second is presented in ^. Section ^ looks at an example of the operation of the 
second algorithm and Section |^ considers some issues of correctness and conver- 
gence relating to the second algorithm, and Section |§| looks at the time and space 
complexity. We conclude in Section ^. 

2. STOCHASTIC AUTOMATA 

In this section we introduce stochastic auto mata using a small e xample. Stochastic 



automata are related to timed automata | Alur and Dill 1994 1, however stochas- 
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tic clock settings are used, instead of the strictly deterministic timings used in 
timed automata. We begin with the formal definition of stochastic automata, then 
present a simple example. We use the definition of stochastic automata presented 



in [D'Argcnio et al. 1998 



Definition 1 . A stochastic automaton is a structure (S, Sq , C, A, — >, k, F) where: 

• >S is a set of locations with sq G S being the initial location, C is the set of all 
clocks, and A is a set of actions. 

• — E> C5x(Ax 'Pfln(C)) x S is the set of edges. If s and s' are states, a is an action 

a,C 

and C is a subset of C, then we denote the edge (s, a, C, s') G — e> by s — e> s' 

a 

and we say that C is the trigger set of action a. We use s — 1> s as a shorthand 

a,C 

notation for 3C.s — E> s' . 

• k : S — y Vfm{C) is the clock setting function, and indicates which clocks are to 
be set in which states, where Va n (C) is the finite powerset of clocks. 

• F : C —y (1Z — > [0, 1]) assigns to each clock a distribution function such that, for 
any clock x, F(x)(t) = for t < 0; we write F x for F(x) and thus F x (t) states 
the probability that the value selected for the clock x is less than or equal to t. 

Each clock x G C has an associated random variable with distribution F x . To 
facilitate the model checking, we introduce a function £ which associates locations 
with sets of atomic propositions. 

£ : S h-> p(AP) 

where AP is the set of atomic propositions. □ 

It is necessary to impose some limitations on the stochastic automata which 
can be used with the model checking algorithms. In particular, we require that 
each clock distribution function F x must have a positive finite upper bound and 
a non-negative lower bound, and must be continuous between these bounds. The 
finiteness constraints mean that there are certain distribution functions which we 
must approximate. We further assume that clocks are only used on transitions 
emanating from states in which they are set. 

As an example, consider the simple packet producer (which is a component in a 
large multimedia specification) in Figure |^. This is written 

({s ,s 1 ,s 2 },s Q ,{x,y,z}, 

{tryagain, cone, send, fail}, — >, n, {F x , F y , F z }) 

where 

— = {(s Q , tryagain, {x}, s ), (s , cone, {x},si), 
(si, send, {z}, s ), (s ,fail, {y}, s 2 )} 

K (s Q ) = {x,y}, k(si) = {z}, k(s 2 ) = {} 

and the distribution functions for clocks x, y and z are 

F x (t) = 2t-i 2 ,if t G [0,1] 
= 0, if t < 
= 1 , otherwise 

Fy(t) - i 2 ,if t G [0,1] 

= 0, if t < 
= 1, otherwise 
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Fig. 2. The packet producer 

and 

F z (t) = *, if te[o,i] 
= 0, if t < 
= 1, otherwise 

as depicted. The horizontal axis measures time, and the vertical axis measures the 
probability of the clock being set to a value less than that time. 

The packet producer starts in location s , and attempts to establish a connection 
with its medium. Three options are possible at this stage. Either the medium allows 
a connection, the medium tells the packet producer to try again or the medium takes 
too long and the connection fails (is timed out). These options are modelled in the 
automaton by setting clocks x and y according to the functions F x and F y . If clock 
x expires first then there is a nondeterministic choice between the transition labelled 
cone (which moves the automaton to state Si) and the transition labelled try again 
(which moves the automaton back to state So). This choice is nondeterministic 
because in reality it would depend on the medium, which we have not specified 
here. If clock y expires first, then action fail is triggered (we say that {y} is the 
trigger set of fail) and the automaton moves to state S2- This corresponds to the 
medium taking too long to respond, and nothing further happens. 

This example has been chosen because it is small enough that we can show, in 
their entirety the set of configurations that our model checking algorithms enumer- 
ate. Thus it can be used to illustrate our two algorithms. But, in addition, we have 
chosen it because it is canonical in the sense that it illustrates the key concepts of 
stochastic automata, e.g. simultaneous enabling of multiple transitions generating 
non-determinism. The reader should also notice that this is a good example of a 
situation in which steady-state analysis is not interesting. Specifically, in the steady 
state, all the probability mass will be in state S2- Thus, the sort of questions we 
wish to ask about such a system are about its transient behaviour, e.g. what is 
the probability of reaching state S2 within a particular period of time and indeed 
this is exactly the type of question we will be able to formulate with the logic we 
introduce in the next section and answer with our model checking algorithms. 
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3. A PROBABILISTIC REAL-TIME TEMPORAL LOGIC 
3.1 The Logic 

In this section, we introduce a simple probabilistic temporal logic. The purpose of 
the logic is to express properties that we wish to check the stochastic automaton 
against. The logic we define allows us to check a range of such properties. 

Recall that the region tree contains nondeterminism, and so we resolve this us- 



ing the notion of adversaries (see for example Baier and Kwiatkowska 1998 1 ) . An 
adversary of a stochastic automaton can be thought of as a scheduler, which re- 
solves any nondeterministic choices which the stochastic automaton must make. 
An adversary may vary it's behaviour according to the previous behaviour of the 
automaton, or it may prescribe that for all non -deterministic choices a particular 



branch is always preferred. See [ D'Argcnio 199Sf] for examples of adversaries. 



We assume that when we wish to model check a property against an automaton, 
we are also given an adversary to resolve the nondeterminism within the automaton. 
(Without this adversary, enumerative analysis would not be possible; the provision 
of an adversary is a prerequisite of model checking.) We can now, for example, 
answer such questions as "Given a stochastic automaton and an adversary, is the 
probability of a send event occurring within 5 time units greater than 0.8?" . 

The syntax of our logic is 

ij> ::= tt | ap | ->V I ^i A -02 | [0i U^ c 2 ] ^ V 

::= tt | ap | -.0 | 0i A 2 

where [0i U^ c fa] — P is a path formula. The path formulae can only be used 
at the top level — they cannot be nested. This is because the model checking 
algorithm we give can only evaluate path formulae from the initial state and is a 
necessary restriction of the current approach. Further: c € N (natural numbers), 
ap is an atomic proposition, p £ [0, 1] is a probability value and ~, ~G {<,>,<, >}• 

We can define a number of derived operators. For example, other propositional 
operators are defined in the usual way:- 

ff = -.tt 

01 V 02 = -'(01 A 02 ) 
01 =>■ 02 = ~>01 V 02 

and we can define a number of abbreviations of a number of temporal operators. 



[O^ c 0] ~ p = 


[tt Ur^C 0] ~ P 


[□^ c 0] ~ p EE 


[-iO^ c -i0] ~ p 


[□0] ~ p = 


[□> O 0] ^P 


[O0] ~ p = 


[C^ O 0] ^ V 


V[01 U^ c <j} 2 } EE 


[01 W^ c 2 ] = 1 


3[0l U^c4> 2 } EE 


[0i u^ c 4> 2 ] > o 


VD0 EE 


V[D0] 


3D0 EE 


3[D0] 


VO0 EE 


V[O0] 


3O0 EE 


3[O0] 



where V and 3 are the branching time temporal logic operators, for all and 
exist 



Emerson 1990]. See [Baier et al. 1999] for similar definitions 



With this syntax, an example of a valid formula that we can check would be 
[tt W<io send] > 0.8 which says that the probability of reaching a send event within 
10 time units is greater than 0.8. 
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3.2 Model Checking 

It should be clear that since we do not allow temporal formulae to be nested we can 
use the following recipe in order to model check a formula ip of our logic against a 
stochastic automaton A. 

1. For each until subformula (i.e. of the form [<p\ IA ~ c 4>2\ — p) in ip perform an 
individual model check to ascertain whether 

A \= [</»i U^ c cj) 2 ] 

2. Replace each until formula in ip by tt if its corresponding model check was 
successful, or ff otherwise. 

3. Replace each atomic proposition in ip by tt or ff depending upon its value in the 
initial location of A. 

4. ip is a now ground term, i.e. truth values combined by a propositional connec- 
tive (-i and A). Thus, it can simply be evaluated to yield a truth value. The 
automaton is a model of ip if this evaluation yields tt, and is not otherwise. 

This recipe employs standard techniques apart from the individual checking that 
A \= [<pi U^c4''A — P an( i tlh s is what our two algorithms address. 

4. THE REGION-TREE ALGORITHM 

In this section we introduce the first algorithm. 

In model checking, we take a temporal logic predicate and seek to establish 
whether it is true for our particular specification. For example, we might try to 
establish whether the above stochastic automaton has the following property: Is 
the probability that a packet will be successfully sent within ten time units greater 
than eighty percent? In order to do this, we need to define a means by which we 
can check the stochastic automaton against this logic. To achieve this the temporal 



logic and the specification must have the same semantic model. In [D'Argenio et al. 



1998], stochastic automata are given a semantics in terms of probabilistic transition 
systems, and so the temporal logic is given a semantics in terms of probabilistic 
transition systems as well, see Appendix |A| 

4.1 Region Trees 

For practical purposes, however, we cannot construct the probabilistic transition 
system, since it is an infinite structure, (both in branching and depth.) We instead 
construct a region tree from the specification. This is finitely branching, but may 
be infinite in depth. Thus, a particular region tree represents an unfolding of the 
stochastic automaton to a certain depth. In fact, we use the temporal logic formula 
to construct a probabilistic region tree, which is used to verify the temporal logic 
formula. More precisely, the region tree is expanded until sufficient probability has 
been accumulated to ascertain the truth or falsity of the formula (this will become 
more clear shortly.) In this section, we describe how to construct region trees from 
stochastic automata. 

We begin with the definition of a valuation, which we use to record the values of 
all the clocks in a particular state at a particular moment in time. The unique clock 
a € C, which we add to the set of clocks, is used to facilitate the model checking. It 
keeps track of the total time elapsed in the execution of the stochastic automaton, 
but plays no part in the behaviour of the automaton. 

Definition 2. A valuation is a function v : C[j{a} — > 7?.1J{-L} such thatv(x) = 
_L or v(x) < x max , where x ma x is the maximum value to which clock x can be set. 
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If d £ 1Z>o, v — d is defined by \/x 6 ClJ{a}.(t; — d){x) 
min(w) returns the value of the smallest defined clock. 



dcf 



v(x) — d. The function 
□ 



Since we assume that clocks are only used in the states in which they are set, 
there is no need to remember their value once the state has been exited. Only 
the clock a maintains its value; the rest are set to J_. At the initialisation of a 
stochastic automaton, clock a is set to some natural number, (we will show later 
how we choose this; it depends on the formula we are interested in) and all other 
clocks are undefined. We define this initial valuation as O n , if O(a) = n. 

We also need a notion of equivalence between the valuations, which will enable 
us to construct the regions within the probabilistic region tree. The issue here is 
the following. Although the size of the tree will be potentially infinite, at each node 
we wish to have finite branching. We achieve this because, although there are an 
infinite number of valuations possible for any particular state, there are a finite 
number of valuation equivalence classes. This gives us the finite branching. 

Definition 3. Two clock valuations v and v' are equivalent (denoted v = v' ) 
provided the following conditions hold: 

— For each clock x £ C[J{a}, either both v(x) and v'(x) are defined, or v(x) =_L 
and v'(x) =±. 

— For every (defined) pair of clocks x,y € C[J{a}.v(x) < v(y) v'(x) < v'(y). 

The same clocks are defined in each valuation, and the order of the values of the 
defined clocks is all that is important. □ 

The reason that the order of the values of the defined clocks is all that is important 
in the definition of a valuation equivalence class is that the actions are triggered 
by the first clock to expire. Therefore we only need to know whether one clock 
is greater than or less than another. Also note that there is a probability of zero 
that different clocks are set to the same value. This is because all distributions are 
assumed to be continuous. 

We are now in a position to describe how a region tree is constructed from a 
stochastic automaton. Intuitively, we build the region tree by "unfolding" the 
stochastic automaton. At each newly reached state, we calculate all possible 
valuations (up to ~) and the probabilities of each one, then from each of these 
(state, valuation) pairs we calculate the possible new states and repeat. 

Suppose we wish to construct the region tree for the stochastic automaton in 



The resulting region tree (up to a particular level of unfolding) is given in Figure ^. 
The first node is labelled with the location Sq, where the SA starts, the valuation 
Oi, (i.e. (1,_L,_L)) since clocks x and y have not yet been set, and clock a is set 
to value one. Clock a is set according to the time value on th e for mula in which 
we are interested; we will give the example formula in Section |4.2| . The clocks x 
and y are then set, giving a potential 3! — 6 different equivalence classes. However, 
these can be reduced to two by observing that clock a will be fixed on 1 and 
Xmax = Umax = 1 and the probability of either x or y being set to exactly 1 is 
zero[|. Using the convention that we subscript the clock variables by the iteration 
number, in order to distinguish different settings of the same clock, the two possible 
equivalence classes are therefore vo(y) < vq(x) < vo(a) and vq(x) < vo(y) < vo(a), 
where vo(a) — 1 in both cases. For convenience, we will write Xq for vq(x), yo for 
vo(y) and an for vo(a). 

1 This coincidence of a, x max and Umax is assumed in order to simplify our presentation; the next 
iteration illustrates the general case. 



Figure 
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Fig. 3. The region tree 



If clock x is set to less than clock y, the automaton will allow time to pass in 
location sq, and each clock will count down, until clock x reaches zero. Then, 
either action tryagain or action cone will fire (the choice is nondeterministic) , and 
the automaton will enter location sq or s\ respectively. The time at which this 
occurs will obviously vary according to the initial value of the clock x. The possible 
locations entered are depicted by regions 3 and 4 in the region tree in Figure ||, 
where clocks x and y (since they are irrelevant in these regions) are not recorded. 
The initial value of clock a when moving from region 1 to either region 3 or region 

4 will be 1 — xo (we will denote this value as ai). Thus, it will be in the range (0, 1). 
If clock y is set to less than clock x (represented by region 2), then the action 

fail fires, causing the automaton to enter location S2, and this is depicted by region 

5 in the region tree. Again, all we can say about the value of clock a at this stage 
is that it lies in the range (0, 1). 

From region 3 there are two possibilities. Either clock z is set to less than a\, 
(region 6), or it is set to greater than a\ (region 7). From region 6 the action send 
will occur before the clock a expires, moving the automaton to location sn and the 
region tree to region 14. From region 7 the clock a will expire before the action 
send occurs. The region tree moves to region 15, and the automaton remains in 
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state s±. 

From region 4 (location Sn) both clocks x and y are reset according to their 
probability density functions, to values X\ and y\. Since we cannot now be sure 
about the value of clock a, we have 3! = 6 equivalence classes, and these are 
represented by regions 8 to 13 when we unfold the SA another level. 

In regions 8 and 9 a\ is less than the (new) initial values of clocks x and y: these 
regions represent the case where clock a expires before either of x\ and y±. When 
we consider a particular temporal logic formula this will represent the case where 
time has run out, and so the region tree moves to either region 16 (if y expired 
first) or region 17 (if clock x expired first). 

Regions 10 and 11 represent the valuation equivalence classes where x\ is less 
than both a\ and y\, and so from these clock x will expire first, either action 
tryagain or cone will be performed, and the stochastic automaton will enter either 
location so or si (regions 18 — 21). 

Region 12 and represent the valuation equivalence classes where y\ is less than 
both ai and xi, so clock y will expire first, action fail will fire, and the automaton 
will enter location S2 (region 22 — 23). 

The region tree can be expanded further if necessary. There is no need to con- 
tinue to expand regions 5, 15, 16, 17, 22 and 23, because in all of these either the 
clock a has expired or the stochastic automata has reached location S2, which is 
a deadlocked state, and there is no further information to be gained. In Figure [|, 
further regions are derived from region 14 in the same way as above; these are 
needed when we build the probabilistic region tree in the next section. 

4.2 Probabilistic Region Trees 

Given a stochastic automaton, adversary and formula ip — [4>i U^ c <j)2\ — P the 
model checking algorithm consists of a number of iterations which are repeated 
until the formula is found to be either true or false. 

An iteration unfolds the region tree by expanding each leaf node. At each itera- 
tion stage there are two steps. The first step resolves the nondeterministic choices 
in the newly expanded region tree using the given adversary. The second step then 
calculates the probabilities on each node in the newly expanded part of the tree. 

The region tree (Figure ||) represents an unfolding of the stochastic automaton 
without the nondeterministic choices being resolved. The probabilistic region tree 
(Figure ^) records the resolution of the nondeterministic choices and the probabil- 
ities at the final nodes represent the chances of taking the particular sequence of 
actions that end in that node. 

At each iteration, we update the information we have on the probability of a path 
satisfying the formula. To do this, we define three new propositions, and each node 
of the probabilistic region tree is labelled with p, f or u: p, if it has passed (it is the 
end of a path which models the bounded until formula ip); f, if it has failed (it is 
the end of a path which cannot model ip), or u, if it is undecided. We also have two 
global variables, Sp and Sf, which keep running totals of the probabilities of the 
pass and fail paths. 

The basic idea of the model checking algorithm is that we check the values of Sp 
and Ef at each stage, and if we cannot deduce from these the truth or falsity of 
the formula we are checking, we look more closely at the undecided nodes. That is, 
we extend the undecided paths by each possible subsequent action, label these new 
nodes p, f or u, and calculate their probabilities. We then add these probabilities 
to Sp and Sf and repeat. 

We will begin by demonstrating the technique for an example. The full algorithm 
appears as appendix Consider the example stochastic automaton (Figure ^). 
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Fig. 4. The probabilistic region tree 



Let us consider the formula 

^ = p VWW<! 2 ]>O.9 

where 0o (resp. </>i, 2 ) is the proposition that we are in state so (resp. s±, S2). 
The question^ we are therefore asking is: is the probability of reaching location s 2 
(failing) within one time unit greater than 0.9? 

Note that a steady state analysis will tell us only that the automaton will fail 
(reach state s 2 ) eventually, but here we want to obtain information about the 
transient behaviour of the automaton. The nondeterministic choice that has to be 
made is between location s± and s 2 . We will consider the benevolent adversary, i.e. 
the one that always chooses location si. 

Consider region 1 first (Figure ||). It has two possible outgoing transitions, and 
the choice between them is made nondeterministically. So we must refer to the 
adversary, which chooses location si, that is, region 3. Region 4 is not generated. 
We note that the value of clock a is greater than zero (so time has not run out), 



2 In fact, the algorithm can easily be adapted to handle questions such as "what is the probability 
(to within some e) of a formula such as [</>o £V<i <f>2] being true?". 
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and that proposition p$ V p\ is true (so the temporal logic formula is able to be 
satisfied), so this region is labelled with u (undecided). 

In region 5 proposition p2 is true, and clock a is still greater than zero, so this 
region is labelled as passed p, and region 5 becomes a terminal node. 

In region 6 a\ is greater than the (new) initial value of clock z, and therefore the 
send action will fire before the clock a expires. The region is therefore labelled u. 

In region 7 a\ is less than the (new) initial value of clock z, and therefore time 
will run out before the send action has a chance to fire. The region is therefore 
labelled f. 

From region 6 the send action moves the automaton to location Sq (region 14), 
and from here there are 6 possibilities for the setting of the clocks. 

Regions 24 and 25 represent the valuation equivalence classes where a\ is less 
than xi and y\. Since clock a will expire before either clock x or clock y, we know 
that these paths will not reach location S2 in less than one time unit, so regions 
30 and 31 will be labelled f. The remainder of the tree is generated in a similar 
manner. 

Figure || represents two unfoldings. In order to determine whether the formula 
is true we also have to calculate the probabilities on the nodes. If the sum of the 
pass and the sum of the fail nodes is sufficient to tell us whether the formula is true 
then we can stop here, otherwise we unfold the tree another level. 

To determine the probabilities on the arcs, we need to use probability density 
functions P x , P y and P z of the functions F x , F y and F z , which we find by differen- 
tiating F x , F y and F z between their upper and lower bounds and setting to zero 
everywhere else. 

P x (t) =2 — 2*, if t e [0,1] 
0, otherwise 

P y (t) = 2t,if t e [0,1] 
0, otherwise 

P,(i) = l,if t e [0,1] 

0, otherwise 

Evaluating the function F x at a point a gives the probability that clock x is set 
to a value less than a, and if a > b, then F x (a) — F x (b) gives the probability that 
clock x is set to a value between a and b, provided a and b are constants. The same 
calculation using the corresponding probability density function (pdf) would be 
J, P x (x)dx, which at first sight appears more complicated. The advantage is that 
these functions can be used to calculate the probability that clock x is set to a value 
less than y, where y is a random variable set according to the distribution function 
F y . If, for example, we wished to calculate the probability of the equivalence class 
in region 1 (vq(x) < v (y) < v (a), where vq(o) = 1) we would evaluate J Q y P x (x)dx, 
to give us a function that returns the probability that v (x) is between and y, 
multiply this by the pdf P y (y), and integrate between zero and one: 

/ / P x (x)dxP y (y)dy 
Jo Jo 

which gives us the probability that x will be less than y, where x and y are random 
variables conforming to the distribution functions F x and F y . 

We will now evaluate the probabilities of some of the arcs in the example. In the 
following, we will continue to subscript the clock variables by the iteration number, 
in order to distinguish different settings of the same clock. 
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1 rl-zi pi 

/ / P y {yo)dyoPx(xo)dxoP z (zi)dzi 

JO Jx 



1 

Py (yo)dyo Px (xo)dx P z (z 1 )dz 1 



Table 1. The integrals 

In our example, to determine the probability on arc (0,2), where the value to 
which clock y is initially set (which we will refer to as yo) is less than the value to 
which clock x was initially set (xq), (yo < xq) we perform the double integration 

1 rX 

/ 2y Q dy Q (2 - 2x )dx 
Jo Jo 

which evaluates to i. 

Arc (0, 1) must have the value 1 — | = |, since it is the only other possibility, 
and can be calculated as 

el 

2y dy (2 - 2x )dx 

These two arcs represent the setting of the clocks, and are therefore instantaneous. 

From Region 2 the only region which can be reached is the leaf node region 5, 
and therefore the arc (2, 5) has probability 1. 

Calculating probabilities on the paths through region 3 is more complicated. Con- 
sider arc (3, 6) first. In fact, we must calculate the probability of the path (0, 1,3,6) 
in its entirety rather than determine separately the conditional probability of arc 
(3, 6). We do this as follows. 

The clock setting information we know is: the first time the clocks x and y are 
set, the initial value of x is less than the initial value of y (xo < yo)] and when z\ 
is set, the sum of xq and z\ is less than the initial value of clock a (xq + Z\ < 1). 
These constraints are captured as the combination of the integrals J x P y (yo)dyo 
(to ensure that Xq < yo < 1), jo * % Px( x o)dxo (to ensure that Xo + Z\ < 1), and 
Jq 1 Pz(zi)dzi (since all constraints have been captured in the first two integrals.) 

The combination is given as the first integral in Table [I] and equals | . 

The path (0, 1, 3, 7) differs only in the fact that a\ (= 1 — Xq) is less than z%, 
and can be calculated as the second integral in Table [l] which equals ^ . The only 
difference is that P x (xq) is integrated between 1 — Z\ and 1. 

At this stage in the algorithm, Ep = | and Ef = Since Ef > 1 — 0.9 we can 
deduce that the formula is false, and in this case, there is no need to unfold further 
the node labelled u. 

The accuracy with which we know the values of Ep and Eu will increase as the 
probabilistic region tree is extended, and in some cases it may need to be extended 
to infinity for perfect accuracy. However, we can achieve accuracy to within an 
arbitrary tolerance e with a finite probabilistic region tree. 

The major drawback of this algorithm is its complexity: with every new unfolding 
of the probabilistic region tree not only does the number of nodes to be considered 
increase, but also the number of integrations required to determine the probability 
on a single node increases exponentially. It therefore becomes intractable after a 
few iterations. This is the issue we try to tackle with the second algorithm. Rather 
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than integrate the probability density functions, we discretise the ranges of the 
functions and work with the resulting approximations. 

5. THE MATRIX ALGORITHM 

In this section we present an overview of the second algorithm. The second algo- 
rithm takes a stochastic automaton SA, together with a bounded until temporal 
logic formula TL, a time step parameter S and an adversary pick. For convenience 
we will present only the case where TL is of the form [4>q U< c <j>i] > P- Minor mod- 
ifications to the algorithm would allow any of > p, < p or < p. We use the atomic 
propositions </>o and <j)\ as part of the formula because anything more complex can 
be reduced to these by standard model checking techniques. Using < c guaran tees 
that the algorithm will terminate, although we discuss the > c case in Section 6.1 . 

A single iteration of the algorithm will return one of three results: true, false or 
undecided. If it returns true, then the automaton models the formula. If it returns 
false, then the automaton does not model the formula. If it returns undecided, then 
the algorithm was unable to determine whether the automaton models the formula. 
In this case, the algorithm can be re-applied with a smaller value for the timestep 

6. The question of convergence to the correct answer as S tends to zero is discussed 
in section [?]. For the remainder of this section we assume S to be fixed. 

A stochastic automaton has a finite number of clocks each with a probability 
distribution function (pdf). For each state, the set of clocks has an (arbitrary) 
order, and the algorithm makes use of this ordering^. We assume that each clock 
has non-zero lower and upper bounds on the values to which it can be set. The 
first of these is a new constraint and was not required for the first algorithm. This 
has been done so that 8 can be initially chosen to be less than the minimum of all 
these lower bounds. 

The algorithm works by creating a snapshot of the automaton at each time point 
7i(5 (n e N)Q and extracting some global information about the probability of the 
formula [4>o U< c (f>i\ being satisfied at this point.^J To build the next snapshot, the 
algorithm picks out at each time point nS the transitions that the automaton is 
capable of during the next interval of length S. Because S is less than the minimum 
of all the clock lower bounds, a maximum of one transition per path can occur in 
each interval. Recording all possible states of the automaton at each time point is 
therefore enough to record all the possible transitions. 

The algorithm stops when either enough information has been gathered to deter- 
mine the truth or falsity of the formula, or enough time has passed so that nd > c, 
and allowing time to pass further will make no difference to the information we 
already have. In this case the result undecided is returned. 

5.1 Data structures 

The principal data structures used by the algorithm are matrices. For each state 
s in the stochastic automaton we derive a matrix for a given time t (which is a 
rational number and calculated as n5), denoted matrix(s,t), which is a record of 
the probabilities of the various combinations of clock values in state s at time t. 

Each matrix matrix{s, t) will have #k(s) dimensions. Each dimension is associ- 
ated with a particular clock, and the ordering of the dimensions corresponds to the 
ordering of the clocks. The dimension associated with a clock c will have [ Cn y" ] 

3 However, the choice of ordering is arbitrary and does not carry any meaning. Any ordering will 
be sufficient. 

4 We will speak of the time instants generated by nS (n S N) as time points. 

5 We also require that 3n.n8 = c, which ensures that one of the snapshots will be at exactly time 
c. 



Stochastic Model-Checking for Multimedia • 15 



entries, where c max is the largest value to which the clock c can be set, and |~ c ™ aj: "| is 
the smallest integer greater than or equal to Cm g ax . For a clock a, we will abbreviate 

The valuation function v gives the value of a particular clock: v(ci) is the value 
of clock 

Each entry in the matrix matrixes, t) is the probability that at time point t, the 
automaton is in state s, and each clock is within a particular time range. Thus, the 
value matrixes, t)[k\ . . . k n ] is the probability that at time point t, the automaton 
is in state s, and v(ci) £ (S(ki — 1), Ski] for each clock c$. 

A further data structure we shall need is live(t), which is the set of states "live" 
at time t (i.e. their matrices at time t contain at least one non-zero entry, and the 
formula is still undecided) . In order to get an accurate picture of the automaton at 
time t + S, we must take into account all states live at time point t. 

A snapshot of the automaton at time t is the set of all matrices matrixes, t) 
where s is in live(t). 

Let pr(ci £ (5(ki — l),Ski]) be the probability that clock a is initially set to a 
value in the range (8(ki — 1), Ski]. Before the algorithm proper begins, we calculate 
all these values from the clock probability distribution functions, which are entered 
into the algorithm as part of the stochastic automaton. 

5.2 Variables 

The algorithm also uses a number of auxiliary variables. 

prob{s, t) is the probability of entering state s during the time range (S(k — 1), Sk] 
(where t = Sk) and is defined for states s live at time S(k — 1), and s' live at time 
Sk. 

new states (s,t) is the set of states which can be reached from a state s during a 
time range (S(k — 1), Sk]. 

totaljpass is a probability value. It is incremented at each iteration. The iter- 
ations of the algorithm correspond to the time points, and totaljpass records the 
probability of the automaton having passed the formula at that time, totals] ail is 
also a probability value; it records the probability of the automaton having failed 
the formula as the algorithm progresses. 

error is an upper bound on the possible errors of totaljpass and totals] ail. 
After an iteration, we know that the actual probability of the automaton having 
passed the formula is in the range [totaljpass, totaljpass + error], and similarly for 
total jf ail. 

5.3 Overview 

The second algorithm is given in detail in Appendix We begin here with a 
pseudocode description. 

build matrix(so, 0) 

check formula against srj and t = — > pass 

-> fail 

| undecided 

repeat 
t:=t + 5 

forall locations s in live(t — 5) 
build matrixes, t) (record possible new locations) 



update live(t) 
forall locations s' in live(t) 



(increment probability of entering new locations) 
( increment error ) 
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check formula against location: 
if pass then add probability to total -pass 
if fail then add probability to total -fail 
if undecided then update matrix(s' ,t) 
until (formula has passed, or 
formula has failed, or 

t has reached the limit set by the formula) 
set all locations undecided at last iteration to false 
if total-pass > formulaprobability then output pass 
elseif total-fail > 1 — formulaprobability then output fail 
else output undecided 



We now present the formula for initially calculating matrices, then describe the 
algorithm in overview, outlining the procedures involved. 

If there are n clocks in state s, then matrix(s, t) is calculated using the probability 
distribution functions of the clocks in state s as follows: 

VI < h < Ni 

n 

VI < k n < N n • matrix(s, t)[ki ...k n ) := J]?w(v(q) <E (<5(fc; - I), Ski]) 



The algorithm begins by calculating matrix(so, 0), where so is the initial state 
of the stochastic automaton. 

live(0) will either be {sq} or the empty set, according to whether the formula 
TL is made true or false by state sq, or whether we cannot yet decide. This is 
determined as follows. If state sq models proposition <j>\, then the formula TL is 
immediately true and live(0) is the empty set. Otherwise, if s$ models (f>o we cannot 
yet decide, and so live(0) contains sq. If the state models neither proposition then 
the formula TL is immediately false, and live(0) is the empty set. 

If the initial step does not determine whether the formula is true or false, we 
perform a number of iterations. Each iteration builds the snapshot at time point 
t + S, based upon the snapshot at time point t. The sequence of snapshots build 
progressively more information as to whether the stochastic automaton has passed 
or failed the formula. 

In the case of a bounded until formula with a < c subscript^, the number of 
iterations is finite (i.e. the algorithm always terminates) because the iterations ter- 
minate either when sufficient information has been extracted to determine whether 
the formula passes or fails, or after the ^th iteration, since the formula cannot 
become true after time c. 

If the information at time t is not enough to determine the truth or falsity of the 
formula, we build the snapshot for time point t + S. We now describe an individual 
iteration. 

An iteration consists of two sections. In the first, we consider all of the states 
which are currently undecided. These are all the states in live(t). For each state we 
create the matrices at time t + S, update live(t + S) and calculate prob(s' , t + S) for 
states s' which can be reached in the interval (t, t + S]. In the second, we look at all 
states which can be reached in the interval (t, t + S], and consider them with respect 



[<j>o I4<c 0i] > P- See Section 6.1 for a discussion of how >c time bounds are handled. 



Stochastic Model-Checking for Multimedia • 17 



to the temporal logic formula. We then either update the global probabilities, if 
the states cause the formula to pass or fail, otherwise we update the respective 
matrices. 

Note that in this algorithm a matrix is updated at most twice. Once within pro- 
cedure newJimejmatrix(veieT to Appendix [c]) , if the state was live at the previous 
time, and once within the procedure new 'state -matrix ', if the state is reachable via 
a transition in the previous interval. 

5.3.1 Creating and updating matrices. We begin with some necessary notation. 
Let us assume S is a fixed rational number greater than zero. 

Definition 4. If c\, . . . ,c„ are the clocks on state s, a valuatior^ is the vector 
of results of the valuation function v{ci) from clocks to 1Z which gives the values of 
each of the n clocks. 

Two valuations v and v' are (5—) equivalent if 

Va.Bkt G N.v(a) G (5(ki - 1), h] A v'(a) G - 1), h] 

A valuation equivalence class (or clock configuration) is a maximal set of equiv- 
alent valuations. □ 

If 6 is understood, we can abbreviate this configuration as (fci, . . . ,k„). For a 
state s and a time point t, the probability Yii=i P r i v ( c i) G (<5(fc; — l),5ki\) is an 
(s, t)-clock configuration probability (or just a clock configuration probability when 
s and t are understood). 

There are two different procedures for updating a matrix. The first (encapsu- 
lated in the procedure new dime -matrix) corresponds to the situation within the 
stochastic automaton where time passes, but the state remains unchanged. In this 
case we must shift the clock configuration probabilities in the previous matrix down 
by one index step (which corresponds to S time passing) and add the result to the 
matrix we are updating. 

We also at this stage determine the new states which can be reached from the 
current state during the S time passing, and the probability of entering these states. 
We do this by looking at all the clock configurations where at least one of the indices 
has the value one. If the clocks are set within such a configuration then we know 
that at least one clock will expire during the ensuing S timestep. 

If only one index in the configuration has the value one then only one clock can 
expire, and only one state can be entered from this clock configuration, and so that 
state is added to the set of states which can be entered from the current state at 
the current time. 

If more than one index in the configuration has the value one, then we simply do 
not go any further into the automaton and the configuration probability is added 
to error. 

The second way to update a matrix corresponds to a transition from one state to 
another within the automaton. It is described in the procedure new state-matrix. 
For each matrix entry we calculate the clock configuration probability, multiply it 
by the probability of moving into this state at this time, and add it to the matrix 
entry we are updating. 

5.3.2 Termination of an iteration. When the iteration terminates, it will out- 
put one of three results: true, false or undecided, true means that the automa- 
ton models the temporal formula, i.e. SA \= [4>q U< c (j>i] > p. false means that 



7 We alter the definition of valuation slightly here for the second algorithm. 
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SA Y= [4>o W< c 0i] > P, and undecided means that the algorithm could not accu- 
mulate enough information to decide whether or not the automaton modelled the 
formula. 

The algorithm makes the output decision based on the three global variables 
total-pass, total-fail and error. 

total-pass is a lower bound on the probability that the stochastic automaton 
models the formula, and total-fail is a lower bound on the probability that the 
stochastic automaton does not model the formula, error is the largest amount by 
which total-fail or total-pass may be wrong. In a sense, it records the size of the 
uncertainty introduced by the choice of 8. 

If neither of these situations holds then the errors introduced by the algorithm 
are too large to determine an answer with this value of 6. In this case, we can 
rerun the algorithm with a smaller 8, and in section |?] we show that the sum of 
the errors tends to zero as 8 tends to zero. Note, however, that in the case where 
the probability that SA models [4>o ^<c4'i\ is exactly p, we cannot guarantee that 
there will be a 8 small enough to allow the algorithm to generate a true or a false. 
This is the sort of limitation that has to be accepted when working with generalised 
distributions. 

6. EXAMPLE 

The second algorithm requires slightly more stringent restrictions on the stochastic 
automaton than the first one, because the clock distribution functions must have 
positive lower bounds, (as opposed to the non-negative lower bounds required by 
the first). Therefore in order to illustrate the second algorithm, we will use the 
automaton in Figure pi but alter slightly each of the clock distribution functions, 
by shifting each of them half a time unit to become 

F x {t) = 2t-t 2 ,ii te (|,|] 

= 0, if t < | 
= 1, otherwise 

F y (t) = i 2 ,if te (|,|] 

= 0, if t < \ 
= 1, otherwise 

and 

F z (t) = t, if te(i§] 

= 0, if t < I 
= 1, otherwise 

In this section, we will consider the temporal formula [(ao V ai) U < za2\ > \, 
where Si ^ a^, i € {1, 2, 3}. 

We now illustrate this algorithm by applying it to the example^. We set 8 equal 
to i. 

Sections A, B and C below correspond to the sections A,B and C in the algorithm 
description in Appendix |c[ Within section C, line numbers correspond to the line 
numbers of the algorithm. 



8 The type of situation where this algorithm would do very badly is if one clock has a very small 
lower bound and all the rest have a very high lower bound. This is accentuated if the first clock 
is hardly used. It might even be that the state where the first clock is used is unreachable or has 
a very low probability of being reached. Thus a criterion for the algorithm to work efficiently is 
that all pdf lower bounds are "similar" . 
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Section A. This section initialises all the variables to zero, and calculates all the 
probabilities of clocks falling in the ranges (0, 6], (6, 26] etc. from the probability 
distribution functions entered as part of the stochastic automaton. 

In our example, the probabilities that the clocks x, y and z are in the ranges 
(0, 6], (6, 26} or (26, 36} are given by 

x y z 
{0,6} 
(6,26] | \ \ 

(25, 35] \ | i 

These are easy to obtain from the clock probability distribution functions. In- 
deed, the ease of determining these probabilities is the main benefit of this algorithm 
and contrasts with the intractable manner in which the integrals explode in the first 
algorithm. 

Section B. The initial state sq does not model ax, but it does model the propo- 
sition an, and so the procedure initjnatrix is called. This returns matrix(sQ, 0) 
which is as follows 
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and is easily derivable from the probabilities above. The procedure also sets live(0) 
to {s }. 

If N x is the upper bound of x, and N y is the upper bound of y, there will be 
\N X x |] entries on the x axis, and \N y x |] entries on the y axis, so in this case 
(where iV x = |, JVj, = § and 6 — |), we get a 3 x 3 matrix. 

This matrix tells us e.g. that when the clocks in the initial state are first set, 
the probability of clock x being set within the range (6, 26] and clock y being set 
within the range (2<5, 3<5] is |. That is, for the clock configuration ((6,26], (26,36]), 
the clock configuration probability is |. 

Section C. We now enter the iterative part of the algorithm, where each iteration 
corresponds to increasing the time by one time unit (6) , and the snapshot produced 
at the end of iteration n corresponds to a view of the automaton at time n6. The 
three global probability values^] are all still zero (lines 1-la), so ct (current time) 
becomes 6. Only the state so is live at time zero, so new -time-matrix is called 
(line 6) for matrix(so, 6). This returns a number of parameters: matrix(so, 6), 
newstates(si, 6),prob and error. 

The procedure new -time-matrix will return the matrix(so, 6) as 
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where each clock has advanced one time unit from matrix(so, 0). So, at time 6, 
the probability of clock x being within the range (0, 6] and clock y being within the 
range (6, 26} is | . 



9 These are the probability values that are updated throughout the algorithm: totaljpass, total-fail 
and error. 
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The probability of staying in state so for at least half a time unit is 1; this 
follows from the fact that no clock can be set to less than S (| time unit). Thus 
prob(s ,8) = 1. 

None of the edge values (those with at least one clock in the range (0, S]) of 
the previous time matrix (matrix(so,0)) is non-zero (so there is no possibility of 
any clock reaching zero and causing a transition to fire). The second half of the 
procedure (lines 10-23, which would determine the new states reached from state 
s ) is therefore not executed and the global probability values (totaLpass, total-fail 
and error) are all still zero. newstates(so, 8) will be returned as {}, since no new 
states can be reached at time 8. 

The next step (lines 7-11 of section C) is to calculate the live states at time 8, 
and since remain(so,5) = true (it is possible to remain in state so at time 8) wc 
include so- 

Since there are no states which can be reached from state so in the time interval 
(0, 8], lines 12-22 of section C are not executed. 

All of the global probability values are still zero, (i.e. we don't have enough 
information to decide the truth or falsity of the formula at this stage, lines 1-1 a 
of Section C), and 28 < 2 (we have more time in which to gain more information, 
lines 2-3 of Section C) , so we begin a second iteration. 

On the second iteration of the while loop, ct is set to 28. Only s was live 
at the last iteration (live(S) = {so}), so at line 6 we call new -time .matrix for 
matrix(so, 28). 

This again returns a number of parameters, e.g. matrix(so,2S) becomes 
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where the entry matrix(so, 25) (1 , 1) is taken from the clock configuration (<5, 26], (6, 26} 
in the previous time matrix matrix(s ,8) and thus the probability of staying in state 
so in the interval (6, 25] is |. However this is not the final version of matrix(so, 26), 
because some of the clock configurations lead to transitions which lead back to state 
s . 

All the other clock configurations ((1,1), (1,2) and (2,1)) in matrix(so,6) lead 
to transitions. Lines 10-22 of procedure new dime -matrix are executed for each of 
these three configurations. 

For clock configuration (1, 1), clock x is (arbitrarily) chosen to fire, and we assume 
that the adversary pick chooses the action cone, leading to state s\. Line 13a of 
the procedure adds state s\ to newstates(so, 26), and prob(s\, 28) becomes | (line 
14). Clock configuration (1, 1) is one where some error may be introduced into the 
algorithm result. Choosing clock x and action cone meant that we go to a state 
where the formula TL can still be true, but choosing the other clock may not lead 
to such a state. We therefore allow for the possible error introduced here by adding 
the clock configuration probability to error, which becomes |. Clock configurations 
(1,2) and (2,1) are dealt with similarly, but error remains constant. 

Now, the new dime-matrix procedure is finished, and lines 7-11 of Section C de- 
termine the value of live{28) which is {so, s\, S2}, because at time 28 the automaton 
may be in any state. 

Lines 12-22 of Section C consider each new state that can be reached in time 
interval (8,28]. State so still allows the temporal logic formula to be true, and so 
procedure new state-matrix is called (line 17). However, prob(so,28) = 0, and 
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therefore matrix(so 1 2S) is not altered. 

State si still allows the temporal logic formula to become true (line 13) and so 
procedure new state -matrix is called (line 17). The probability of entering state 
s\ in this interval is |, so matrix(si,2S) is 





3 3 
U 8 8 




12 3 2 



In state S2 the formula is true, and so prob(s2, 25) (g) is added to totaljpass (line 
14). 

In the final iteration, the global probability values become: total-pass = g, 
total-fail = | and error — |. The iterations stopped because the value of time 
became too large — not because the global probabilities contained enough infor- 
mation to make a decision. This means that totaljpass (g) is a maximum possible 
probability value of the formula [(an V a±) U < 3 a{\ (with any clock ordering) and 
totaljpass — error (— g) is a minimum possible probability value. 

Thus, since we wish to determine whether the actual probability value is greater 
than i, the algorithm will output fail. 

If we were interested in a similar formula with a probability value in the range 
[0, g], we could reduce the size of 5, and take snapshots (e.g.) every j time unit. 
This (for the reasons outlined in Section will reduce the size of the error variable. 

6.1 Unbounded until formulae 

As just presented the second algorithm only handles until formulae of the form 

[4>i u< c 4>2\ — p 

however a combination of the second and first algorithms yields a method to verify 
unbounded until formulae, i.e. those of the form 

[01 U >c (j)2] ^ P 

The basic idea is to observe that verification of a formula such as <j)\ U >c 4>2 can 
be split into a conjunction of separate verifications 

(a) Check that 4>i holds at all times until c time units have elapsed; and 

(b) Check that there exists an X > c such that 4>2 holds at time X, and that for 
all times strictly greater than c and less than X, <f>i holds. 

Thus, we can model check formulae such as [cf>i U >c <j>2\ — p in the following way. 

(i) Run (the obvious slight adaption of) the second algorithm to check that (a) 
holds. This will finish with a certain amount of probability mass in the variable 
totaLfail and no probability mass in totaLpass. The reason for the latter is that 
pass states can only be revealed once time has passed beyond c. In addition, 
live(c) will indicate the locations that are still undecided, i.e. from which we 
must explore further. 

(ii) Run the first algorithm using live(c) as the starting locations and the initial 
timing regions determined from the remaining matrices (this can be done in a 
straightforward manner). However, notice that running the first algorithm in 
this situation does not incur the problems of intractability that it does in the 
general case. Specifically, since the time bound on the until has been satisfied we 
ostensibly only have an untimed until verification. Consequently probabilities can 
be assigned to nodes without requiring the global clock to be taken into account 
and thus, they can be evaluated "locally" . Hence, the exponential explosion in 
the number of integrals to be considered does not occur. 
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7. CORRECTNESS AND CONVERGENCE 

For a single run with fixed S, we wish to prove two things: that the algorithm 
terminating with pass implies that the automaton models the formula, and that 
the algorithm terminating with fail implies that the automaton does not model the 
formula. 

If the algorithm outputs pass then the variable totaLpass must be greater than 
p (where p is taken from the temporal formula [<po U< c 4>i] > p). The only place 
where totaLpass gets incremented is line 14 of section C (see full algorithm in 
Appendix ^|). If the current state q models (f>i (and all previous states in the path 
model <po) we add the probability of entering the state q at the current time point. 
If the sum of these probabilities is greater than p then the algorithm outputs pass. 

We will consider the case when the algorithm outputs pass. Consider the initial 
state. Note that for any clock configuration, the probability of all paths which 
commence with the clocks being set somewhere within this configuration is equal to 
the clock configuration probability. Furthermore, for an arbitrary state s and time 
c and configuration, the probability of all paths which go through this configuration 
at this time is the probability of the configuration multiplied by the probability of 
reaching that state at that time. 

The probability of reaching state s at time c is the second parameter passed to 
the procedure newstate_matri^\ 

If every valuation in a configuration corresponds to the same automaton transi- 
tion, and this transition is the final one in a path which models the formula, then we 
add the clock configuration probability (multiplied by the probability of reaching 
that state at that time) to totaLpass. 

This is the only way in which the algorithm adds to the variable totaLpass. Since 
the algorithm only outputs pass if totaLpass is greater than the formula probability 
p, it is clear that the algorithm will only output pass if the automaton models the 
formula. 

If more than one clock in the configuration is in the range (0, 5] then more 
than one of the clocks will have reached time in the interval we are considering, 
and so the clock configuration probability is added to error (line 12 of procedure 
new_timejmatrix) . 

A similar argument applies in the case where the algorithm outputs fail. 

Therefore the algorithm is sound in the sense that if we are given a definitive 
answer, this answer is correct. There remains, of course, the question of convergence 
to the correct answer, and the following theorem summarises the situation. 

Theorem 1 . For every automaton SA and propositions <j)Q and <pi it is the case 
that if SA models [<fio U < c cfti] with probability p, then for any error e greater than 
zero, there is a timestep S greater than zero such that for the formula [4>o U < c 4>{\ > 
q, the algorithm will only return undecided if q G [p — e,p+e\. 

First note that n independent single variable continuous probability distribution 
functions f\...f n can always be combined to give a single n variable probability 
distribution function which is continuous in all dimensions: f(x\ . . . x n ) — fi(x\) x 

' ' ' x fn(.Xn). 

For convenience, consider a location with two outgoing transitions and two clocks 
x and y with distribution functions f x and f y . Because f x and f y are both contin- 



1[ 'ln fact, it is greater than or equal to this sum, because some routes through the transition system 
may have passed or failed the formula already, and therefore would be considered no further by 
the algorithm. 
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Fig. 5. Upper bound on error with clocks x and y. 

uous, if we set f(x,y) = f x (x) x f y (y) we can (by the note above) say that 

Ve > 036 > 0.f(x, x + 5)- f(x, x - 8) < e 

We will show that for any desired size of error we can choose a suitably small 
timestep. 

Now, J™ f(x, x + 5) — f(x, x — S)dx [^] (the probability of the clock valuation 
falling between the two 45 degree lines in Figure ||) is greater than the sum of 
all contributions to the error variables (represented by the squares in the figure). 
Since the number of locations in the stochastic automaton is finite (say N s ) and (for 
bounded until formulas with less than subscripts) the maximum number of visits to 
any location is finite (say N v ) for any desired error e we must ensure that, for every 
location, for the multivariate function associated with that location, we choose e 
such that e < N ^ N . If the timestep is set to the smallest S necessary to ensure 
that every location provides errors less than N ; then total error provided by 
one location (over all time) will be less than jj- and the total error provided by all 
locations will be less than e. 

8. COMPLEXITY MEASURES 
8.1 Time complexity 

The time complexity of the algorithm discussed in Section ^ depends on a number 
of factors, namely 5, t, m, «2 and | S \. The explanation of these parameters is as 
follows: 

— t is the value of time given in the time-bounded until formula: [a U<tb] ~ p] 
— 5 is the chosen timestep; 

-| S | is the number of states in the automaton; 
— rii is the largest number of clocks in a single state and 
— U2 is the largest (positive finite) upper bound of all the clocks. 

An upper bound on the number of matrices which need to be built in a single 
iteration is | S \ , where S is the set of all states in the automaton. 

To calculate the time complexity we also need to calculate the size of the largest 
matrix. Each matrix is multi-dimensional, and 5j£ will be the maximum number 
of entries over all matrices and all dimensions. For example, in the example in 



1 m = min{x max, Umax}, where x max is the largest value to which clock x can be set. 



24 • J. Bryans, H. Bowman and J. Derrick 

Section ^| all the matrices had 2 dimensions and the maximum number of entries in 
any dimension was 3 since 5 = | and ri2 = §. 

An upper bound on the size of the largest matrix will therefore be the number 
of elements in the largest dimension, raised to the power of the largest number of 
clocks on a single state. 

The time complexity is thus bounded by the time taken to update all the possible 
matrices in each iteration of the while loop in the algorithm, multiplied by the 
maximum number of iterations the algorithm will perform in the worst case. This 
latter value is |, therefore the time complexity is 

|x(^rx|5| 

Although this is exponential, the exponent ri\ is something which should in gen- 
eral be fairly small (< 3) because we only allow clocks to be used from the state in 
which they are set. 

In fact, the algorithm could be optimised to provide a better time complexity, 
by limiting the size of the matrices to min(4, ^f-) since there is no need to consider 
the operation of the clock beyond the limit set by the time bound on the temporal 
formula. The size of the largest matrix would therefore be less than (min(|, ^)) ni , 
where ri\ is the largest number of clocks in a single state. 

An upper bound on the time complexity would therefore be 

JxW^fx \S\ 

The time complexity also relies heavily on 5, and the bigger the S the lower the 
time complexity. To see the relationship with 5, note that the upper bound can be 
rewritten as 

(rj xtx(n 2 rx|5| 

8.2 Space complexity 

An upper bound on the space complexity will be proportional to the product of 
the size of the biggest matrix and the largest number of matrices which need to 
be stored at one time. The size of the largest matrix is less than (-nf-)™ 1 , (from 
time complexity calculations) and the largest number of matrices which need to be 
stored at any one time is twice the number of states in the automaton, 2 x | S \. 
The upper bound on space complexity is therefore 

2x(|)"'x | S | 
9. CONCLUSIONS AND FURTHER WORK 

In this paper we have presented two algorithms for model checking bounded until 
formulae against stochastic automata. Both of these algorithms allow systems to 
be described using continuous probability distributions, and we believe that this 
represents an important advance. 

The principal advantage of the first algorithm is its generality: the clocks may 
be set according to any function, providing the corresponding probability density 
function is integrable. The major drawback of the algorithm is its complexity: with 
every new unfolding of the probabilistic region tree not only does the number of 
nodes to be considered increase, but also the number of integrations required to 
determine the probability on a single node increases exponentially. 
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The principal advantage of the second algorithm is its efficiency: the discreti- 
sation of the probability functions means that the calculations required are con- 
siderably simpler. A limitation in comparison to the first algorithm is that the 
probability distributions must have a finite lower bound. 

In addition, an advantage of both the algorithms is that, since the "complete" 
model is at no point generated, the state space explosion (which typically hinders 
model checking) is contained. In particular, all data structures apart from those 
which reflect undecided nodes (i.e. u labelled regions in the first algorithm and live 
locations in the second algorithm) can be deleted. In this sense the algorithms yield 
a form of on-thc-fly exploration - only keeping information about the "leaves" of 
the exploration tree. 

Further work on the second algorithm will include relaxing the restrictions im- 
posed on the stochastic automata, particularly the ability to set and use clocks 
anywhere in the automaton. Being able to do this would allow parallel composi- 
tion. 

It would also be good to increase the expressiveness of the logic, allowing nested 
untils or "greater than" queries, and to extend the model checking algorithm itself 
to allow queries such as "what is the probability of [<j>Q U< c (pi]T' and receive a 
probability value for an answer. 
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APPENDIX 

A. SEMANTICS 

A.l Probabilistic Transition Systems 

The definition of the semantics of stochastic automata is given in terms of prob- 
abilistic transition systems. The definition of probabilistic transition systems is 



reproduced from [D'Argenio et al. 1995 



N is the set of non-negative integers. R is the set of real numbers, and R>o the 
set of non-negative reals. For n G N, let R" denote the nth cartesian product of 

R. R° d =l f {0}. 

A probability space is a structure (Q,J-,P) where f2 is a sample space, T is a a- 



algebra on Q and P is a probability measure on J-. In this work, as in [D'Argenio 
et al. 1998| , we consider only probability spaces isomorphic to some Borel space 



defined in a real hyperspace, whose coordinates come from independent random 
variables. We denote by 1Z(Fi, . . . F n ) the probability space (R n , BiJEV 1 ), P n ) where 
S(R n ) is the Borel algebra on R n and P n is the probability measure obtained from 



F\ . . . F n , a given family of distribution functions. See [Shiryayev 1984 for details. 



Let V = (f2, T, P) be a probability space. Let V : £1 — > O/ be a bijection. We lift V 

to subsets of fl: V{A) d = {D{a) \ a e A} and define P =' {D{A) \ A £ J 7 }. Now, 

it is clear that V(V) d = (ft' ' ,T' ,P o V^ 1 ) is also a probability space. Since V(V) 
is basically the same probability space as P, we say that T> is a decoration and we 
refer to T>(V) as the decoration ofP according to T). This is used when we come to 
give a semantics to stochastic automata. 

Definition 5. Let P(H) denote the set of probability spaces (Q,J-,P) such that 

il C H. A probabilistic transition system is a structure T = (£, Cn, C, T, >) 

where 

(1) E and £' are two disjoint sets of states, with the initial state (To6S. States 
in £ are called probabilistic states and states in £' are called non- deterministic 
states. 

(2) C is a set of labels. 

(3) T : £ — > f (£') is the probabilistic transition relation. 

(4) — >C £' x C x £ is the labelled (or non-deterministic) transition relation. 
We use a' a to denote (a 1 , 1, a) G — >, a' -/-^ for Sa.a 1 —*—> a and a' — > a 
for 3l.a' -U a. □ 
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Since we are interested in timed systems, we set C = A x R>o, where A is a set of 
action names. A timed action transition will be described as a(d), which indicates 
that the action a occurs exactly d time units after the system has been idling. 

Definition 6. A valuation is a function v : C — > R>olJ{-'-} sucn that v(x) < 
%max, where x max is the maximum value to which clock x can be set. The set of 
all valuations is V. If d G R>o, v — d is defined by Vx e C.(v — d)(x) = v(x) — d. 
We assume the set of clocks is ordered so, if C C C, we can write C for the ordered 
form of C and C(i) for the i-th element. Let C C C, n = #C, and £)S R™. We 
define v[C<-^D] by 



«[C<-4D](:b) = 



45_ f J D (i) if x=C(i):f° r some ie{l,...,n} 
_L otherwise 



□ 

This definition will be used when we explain how clock values change as states 



change. It differs from the definition given in DArgenio et al. 1998] because there 



clocks not in the set C maintain their values through this operation. This is because 



in [DArgcnio ct al. 1998 1 clocks may be used to trigger actions in any state, not 
just the state in which they are set. In this work, however, in order to simplify the 
model checking, we insist that clocks are only used in the states in which they are 
set, and therefore there is no need to remember their value once the state has been 
exited. 

The main obstacle now in constructing the probabilistic transition system seman- 
tics is in showing how the clock probability functions are used to construct the 
probability spaces. We do this by defining a decoration function, discussed in Sec- 
tion 



A.l 



Let SA = (S, sq,C, A, — t>,K,F) be a stochastic automaton. Let s be a location 

in S and n = #k(s). Let v be a valuation in V. Let V = {v [k(s)<-^D] \D& 
R"} C V. We define the decoration function V s v : R" -> {s} x V x {1} by 

V S V (D) = (s,v[k(s)*-]D], 1). Notice that V% is a bijection. In the next definition, 
we use the probability space 1Z(F X1 , . . . , F Xn ) decorated according to some T> s v . 

Definition 7. Let SA = (S, sq,C,A, — >, k, F) be a stochastic automaton. The 

actual behaviour of SA is given by the PTS I(SA) = ((5 x V x {0}), (S x V x 
{1}), (sq, 0, 0), A x R> ,T, — >), where in the initial valuation clock a is set to 
some natural number (chosen according to the PRTL function, see Section^), and 
each other clock is undefined. T and ► are defined as follows: 



k(s) — {xi , x n } 

T(s,v,0)=V%(1l{F xl ,...,F Xn )) 



Prob 



a, {a:} 

s —i> s' A deR> A (v - d)(x) < 
W € [Q.d).W.s b '-i s'.jv - d')(y) > 



(s,u,l)^3(s',(«-d),0) 



□ 



Within a stochastic automaton, two forms of uncertainty may arise. One is the 
probabilistic uncertainty associated with the clock-setting. Although we know 
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which clocks are to be set, the choice of values for these clocks is probabilistic. 
This is where the stochastic element of the model arises, and is defined by rule 
Prob. The other is the nondeterministic uncertainty that arises if two actions 
are simultaneously able to be performed, and is defined using the rule Act. This 
nondeterminism is resolved using an adversary (Definition [l0|) . 

Definition of a PTS-path: 

Definition 8. A PTS-path is a finite or infinite sequence of states 

(a , 00,01,01, . . .) 

where, cto is the initial state, for each a' iy there exists a probability space (S,J-,P) 
such that T(cr l ) = (S, T, P), a' t G S and a[ — ► a i+1 . □ 

Definition 9. An SA-path is a finite or infinite sequence 

((s , wo), (so,v' ), (si,fi), {si,v[), . . . , (s n ,v n ), {s n ,v' n ), . . .) 

such that 

— Vq means no clocks are set. 

—v'i 6 TZ{F X1 ,...,F Xn ) where T( Si , v h 0) = V s v (K(F Xl , . . . , F X J). Each valuation 
vl is a possible result of the clock setting functions. 

— (si, v[, 1) (sj+i, f j+i, 0) for some d. Timed action transitions must be allowed 
by the SA. 

— Finite paths end on a probabilistic state. 

□ 

An SA-path is like a run of the SA expanded with clock values. 

Definition 10 . An adversary of an SA is a function mapping sequences of states 
to states 

adv :< s , si, . . . , s„ > — > s n+ i 
such that < so, si, . . . , s„, s„+i > is a run of the SA. □ 

Note that adversaries do not make any reference to time. 

With an adversary, an SA becomes deterministic. The corresponding PTS contains 
no nondeterminism either. 

If 

o = ((s ,0), (s ,v' ), (si,«i), (siX), ■ • ■ , (sk,v k ), (sk,v' k )) 

is a finite SA-path, then a[i] — Si and o(x) is the state at time x. 

1Z(F Xl , . . . , F Xn ) is the Borel space (R™, Z3(R"), P n ) where P n is the unique proba- 
bility measure obtained from 1Z(F X1 , . . . , F Xn ). 

Now, for all j < k, set Aj to be the maximal set of valuations equivalent to Vj 
which lead to state Sj+%. 

Let 

C(s , A , si, ... , Sfe_i, Ak-i, s fc ) 

denote the cylinder set which contains all paths starting at so and going through 
all states Sj(j < k) and valuation sets Aj(j < k). 
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The probability measure Pr on T{Path{so)^^ is identified by induction on k by 
Pr(C{s )) = 1 and for k > 0: 

Pr(C(s Q , A , ...,A k , Sfc+i) = Pr(C(s , A), • • • , **)) ■ f(^fe) 

where -P(-Afc) is the probability of the set and is taken from the relevant Borel 
space. 



A. 2 PRTL Semantics 

In this section, we introduce the semantics for the temporal logic PRTL. 

To facilitate model checking, we use Probabilistic Transition Systems as a semantic 
model for the definition of PRTL. But in order to do this we must resolve two 
problems. The first is that PRTL is a real-time logic — it enables reference to 



specific instants in time — and the abstract definition of PTSs D'Argenio et al. 



1998] docs not contain reference to time. This is easily solved — we simply use 
the PTS generated by a Stochastic Automaton. This contains much more detailed 
state information, in particular, the values of clocks. 

The second problem is that the PTS contains nondeterministic information, and this 
nondeterminism must be resolved before we can use the PTS to assign a semantics 
to our logic. We do this using adversaries. 

Recall the syntax of PRTL: 

ip ::= tt | ap | -.V I A^ 2 | [(pi U ^ c (p 2 ] ^ P 
<j) ::= tt j ap | -.0 | 0i A (p 2 

where c 6 N, a is an atomic proposition, p £ [0,1] is a probability value and 
{<,>,<,>}■ 

The path formulae ip can only be used at the outermost level — they cannot be 
nested. This is because the model checking algorithms only evaluate path formulae 
from the initial state. 



Definition 11. If SA = (S, s ,C, A, — >, k, F) is a Stochastic Automaton and 
PTS = (E, E', (To, C, T, — ►) is the resulting Probabilistic Transition System, then 
E(= E') C,SxV, £ C Ax R>o and <tq — (sq, 0). We must also introduce a function 
£ which maps SA locations to the logical propositions true in that location. □ 

We only need to use the probabilistic states to define the logic, since once a proba- 
bilistic state has been entered the behaviour of the automaton is completely deter- 
mined until the first clock expires. 

The simple formulae <p are defined in the conventional way for each probabilistic 
region a', but the until formulae ip are defined only for the initial region <7n. The 
model checking algorithm does not yet allow path formulae to be established for an 
arbitrary region. 

• s |= tt 

• s |= a, provided a S £(s) 

• s |= <pi A 02, provided s (= <f>\ and s \= 4>2 

• s \= -i(f), provided s ^= 4> 



Path(so) is all paths possible from so, and J r (Path(so)) is the smallest a— algebra on Path(so). 
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If a is an SA-path, and ip a path formula then 

• cr |= [fa U fa] iff 3fc > 0.(a[k] \= fa A V0 < i < k.a[i] \= fa) 

• a \= [fa U ~ t fa] iff 3a; ~ t.(a(x) \= fa A Vy e [0, x).a{y) \= fa) 

and 

• PTS (= [fa U~ t fa] ^P iff Prob(s ,fa U ~ t fa) ^P where Prob(s ,fa d = Pr{p e 
Path(s ) | p |= V} 

Therefore, the Probabilistic Transition System PTS models the PRTL [fa U ^tfa] — 
p provided Prob(s , fa U ~tfa) — P- 

B. FIRST ALGORITHM 

Here, we give the definition of the first model checking algorithm for bounded until 
formulae. We will consider a PRTL formula of the form [fa U < c fa] > p. "less 
than p" queries may be handled in a similar way. 

Assume an adversary Adv, and that each SA location is mapped to either fa or 
and to either fa or ->fa. Note that the algorithm can easily be extended to 
the more general case where locations contain set of atomic propositions. 

Add the (new) clock a to the set of all clocks. 
Construct the PRG node (sn,O c ). 
Set s = so- 
li s fa then stop with no, else 
REPEAT 

For each possible valuation equivalence class [vi] from n(s)\J{a}, form the node 

(*,M)- 

For each new node (s, [vi]) choose a subsequent non-deterministic node (sj,±) 
according to the adversary Adv. 

For each new non-deterministic node (sj, _L) 
label 'p' if Sj \= fa and v(a) > 0. 
label 'f if Sj \£ fa or Sj y= fa or v(a) < 0. 
label 'u' otherwise 

For each node labelled with either 'p' or 'f, calculate the probability of the cor- 
responding path. 
If T, p pr(s, [v]) > p then stop with yes. 
If £fpr(s, [v]) > 1 — p then stop with no. 
Otherwise, repeat for each node labelled V. 

C. SECOND ALGORITHM 

In this section we present a detailed description of the algorithm. It is divided into 
Section A (which initialises variables), Section B (the initial part of the algorithm) 
and Section C (the iterative part). Procedures used are described at the end. 

The lines of code are prefaced with numbers, and the comments are delimited 
with double stars. 

** Section A** 
Model .check (S A, Formula, S, pick) 

** note that the function pick is the adversary, used in procedure new -time -matrix. 
** We are assuming a TL formula of the form [a U < t ai] > p. ** 
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** The > p could easily be changed; the < t is hardwired into the algorithm. ** 

** ** 

** We begin by initialising variables.** 
** ct: (integer) current.time** 

ct := 

** total jpass and totaLfail are reals in [0, 1]. ** 

** At any point in the algorithm, totaLpass is the accumulated ** 

** probability of all the passed paths and totaLfail is the accumulated ** 

** probability of all the failed paths. We initialise them both to zero.** 

totaLpass := 
totaLfail := 

** error is a real in [0, 1]. It is the accumulated probability of all paths ** 

** which, because of the discretisation of the algorithm, we cannot determine exactly.* 51 

** This is where the revised version of the algorithm differs from the initial one.** 

** It is initialised to zero. ** 

** ** 

error := 

** prob(s,t) is the probability of moving (from anywhere) to location s ** 
** at time t. (i.e. in interval (t — S, t].)** 

** For all combinations of locations and times, we initialise prob ** 
** to zero. ** 

Vs G S.\/i < n. 
prob(s, Si) := 

** remain(s,t) is a boolean which is true if the probability of remaining ** 
** in location s during time interval (t — 6,t] is non-zero, false otherwise.** 
** They are all initialised to false.** 

Vs G S.\fi < n. 

remain^s, Si) := false 
** live(t) is the set of locations "active" at the end of ** 
** interval (t - 6,t], which ** 

** we need for calculating the information for the next time interval. ** 
** For all time values, we initialise live to the emptyset. ** 

Vi < n. 

live(Si) := 

** We initialise all values in all matrices to zero.** 
** The are n s clocks in location s.** 
Vs G S. 

V0 < j < n. 
VI < k < N-l 

VI < i Us < N ns . matrixes, 5j)[i\ . . . i ns ] := 

** call procedure for calculating probabilities of clocks falling in the ranges ** 
** (0,(5], ((5,2(5] etc. This comes directly from the clock PDFs, ** 
** and is only calculated once. It is needed for determining the clock** 
**probabilities. ** 

**C is the set of all clocks and F is the set of clock probability functions** 
** This procedure returns pr, which is needed in newstatejmatrix ** 



** and init-inatrix. ** 
clock_config_probs(C, F, 6, pr) 



32 • J. Bryans, H. Bowman and J. Derrick 



** Section B** 

** Consider initial location of SA: s_0 ** 

** If s_0 |= a_l then formula is trivially true. ** 

if s_0 |= a± then 

totaljpass := 1 
** If s_0 |= a_0 then formula is undecided and we must ** 
** unfold SA further. ** 
elseif s_0 |= an then 

** Build the initial matrix, i.e. matrix(sS), 0). ** 

**This will then contain the probabilities ** 

**of all the different clock settings for location s_0 at time zero. ** 

init-matrix(matrix(s-0, 0)) 

** The only location "live" at time zero will be s_0. ** 
live(0) := {s_0} 

** If s_0 does not model a_0 or a_l then formula is trivially false. ** 

else 

total _fail := 1 
end if 



** Section C** 

** Each iteration of the following loop unfolds the automaton by ** 
** one time step of S. States which cause the formula to ** 
** pass/fail are pruned from the tree, and their probabilities added to ** 
** total -pass /total -fail, while the undecided states are recorded ** 
** for the next iteration. ** 

** We continue while the values of totaljpass, total-fail and error ** 
** are not enough to determine whether the formula is true or false ** 
1: repeat 

** Increment current.time ** 
2: ct := ct + 5 

** for all states s that were live at the last clock tick ** 
4: Vs G live(ct - S) 

** set current_state to s. ** 
5: cs := s 

** The procedure new -time-matrix returns ** 

** matrix(cs,ct): the matrix for the current state at the current time. ** 
** It also ** 

** updates the function prob with the probability of remaining ** 

** in the current state at the current time and the probabilities of ** 

** moving to different states at the current time. ** 

** It also updates the value of error. ** 
6: new-time-matrix{matrix{cs, ct), newstates(cs, ct), remainics, ct),prob, error) 

** If the probability of remaining in current state at current time is zero ** 
7: if remain(cs, ct) — false then 

** current state is not live at current time and ** 

** only the states which can be reached from current state at current time ** 

** are added to those live at current time ** 
8: live(ct) :— live(ct) [J newstates(cs, ct) 

9: else ** remain(cs, ct) = true ** 

** The current state, plus all states which may be reached from it at ** 

** the current time, must be added to the live states. ** 
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10: live(ct) := live(ct) 1J{ CS } U newstates(cs, ct) 

11: end if 

11a: end forall ** Vs e live(ct - S) ** 

** Now, we have live(ct) and prob(cs,ct) for all cs in live(ct) ** 

** i.e. all the states we could be in at time ct, and the probability of ** 

** actually entering them in the previous time interval. ** 

** For every state which can be reached at the current ** 
** time, we must see if it causes the formula to pass or fail, in ** 
** which cases we adjust the values for totaljpass or ** 
** total-fail and remove the state from the live set. If we cannot yet ** 
** tell whether the formula is true or false, we must build the state/time matrix. ** 
12: Vq 6 live(ct) 

** if q \= ai, then formula is true ** 
13: if q \= a\ then 

** totaljpass is incremented by the probability of entering q ** 
** from the current state at the current time ** 
14: totaljpass := totaljpass +prob(q, ct) 

** State q is removed from the live set ** 
15: live(ct) := live(ct) \ {q} 

** Otherwise, if q \= an (and q is not a terminating state) ** 
** then the formula may still be true, ** 

** so we must build matrix(q,ct) and keep state q in the live(ct) set. ** 
16: elseif q \= ao A q terminating states then 

** The procedure new state .matrix returns ** 

** matrix{q,ct): the matrix for state q at current time, and requires ** 
** prob(q, ct): the probability of entering state q from the current ** 
** state at the current time. ** 

17: new state jmatrix(matrix(q, ct) 1 prob{q 1 ct)) 

18: else ** If q does not model a_0 or it is a terminating state and also ** 

** it does not model a_l then the formula is false ** 
** total-fail is incremented by the probability of entering q ** 
** from the current state at the current time ** 

19: total-fail :— total-fail + prob(q, ct) 

** State q is removed from the live set ** 

20: live(ct) := live(ct) \ {q} 

21: end if 

22: end forall ** for all states in live(ct) ** 

23: until totaljpass > p ** formula has passed ** 
24: or 

25: total-fail > 1 — p ** formula has failed ** 

26: or 

27: (error > 1 — p A error > p) ** no possibility of a pass or a fail ** 

28: or 

29: ct = t ** time's up.** 

30: if [ct = t) then 

** All states undecided at the last iteration are now false, so ** 

** total-fail is set to 1 — totaljpass - error ** 
31: total-fail := 1 — totaljpass — error 
32: end if 

**** 

** Output result, based on the values of** 
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** totaLpass, total. fail and error ** 
33: if totaLpass > p then 

** SA models formula ** 
34: output pass 

35: elseif ** total. fail > 1 - p ** 

** SA does not model formula ** 
36: output fail 

37: else ** errors are too large; cannot decide ** 

38: output undecided 
39: end if 



** This procedure builds the initial matrix. ** 

** We assume there are n clocks associated with this state, ** 

** and cf is the Zth clock. ** 

** We abbreviate \ upper. bound(c s l °)~\ .y by TV/. ** 

procedure init.matrix(rnatrix(so, 0)) 
begin procedure 
VI < h < JVi 

n 

VI < i n < N n .matrix{s ,Q)[ii ...i n ] := Y\p r ( c i" e fa _ 
end procedure 



procedure new.time.matrix(matrix(cs , ct), new.states(cs, cCj^remainics, ct),prob, error) 

** This procedure updates a matrix by incrementing time, not by ** 

** changing state. We can do this by considering the values in the previous time ** 

** matrix. It also updates the function prob,** 

** and the variable error.** 

** There are n clocks in state cs.** 
begin procedure 
1: VI < h < Nx 

2: VI < i n < N n . 

** If one of the matrix indices is at its maximum value, then the ** 

** probability value in this position must be zero. This is ** 

** because this procedure is always the first to update a state/time matrix. ** 

3: if 31 < n • i t = N t then 

4: matrix(cs,ct)[ii, . . . ,i n ] := 

** otherwise the values in the matrix can be updated simply from the ** 

** values in the previous time matrix. ** 
5: else ** all clocks q are > 1 and < N t ** 

6: matrix(cs,ct)[ii, . . . ,i n ] := 

7: matrix(cs, ct)[i\, ...,«„]+ matrixes, ct — 5)[ii+i, . . . , i n +i] 

** we record the fact that it is possible to remain in this state ** 
** at this time. ** 

8: remain(cs, ct) :— true 

9: end if 
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9a:end forall 

** We now pick out the positions in the previous time matrix which, ** 
** when moved forward one unit in time, result in a new state. ** 

10: VI < h < Nx 

n-.vi <i n <N n 

** If more than one of the previous time matrix indices is one, we know that ** 
** more than one of the clocks will have reached zero by ct, and so we ** 
** add the probability to error. ** 

11a: if #{q I q = 1} > 1 then 

12: error ' — error + m<xtrix{cs , ct — <5) [ii , . . . , i n ] 

12a: else if #{q | q = 1} = 1 

** Given the stochastic Automaton SA, the state cs and the clock cc ' 
** s' is the resulting state. If the clock is associated with more than ** 
** one transition the function pick (the adversary) chooses the ** 
** resulting state. Otherwise the state is the one determined by the ** 
** transition relation of the SA. ** 

13: s' :— pick(SA, cs, cj) 

13a: new_states(cs, ct) := newstates(cs, ct) U{ s '} 

** the probability of entering s' at time ct ** 
** is incremented by the matrix probability ** 

14: prob(s', ct) := prob(s' , ct) + matrix(cs, ct — S)[i\, . . . , i n ] 

22: end if **line 11** 

23: end forall 

24:end procedure 



** This procedure builds a new matrix, where the state is new rather than the time ** 
** We assume there are n clocks associated with this state, ** 



** and cf is the Zth clock. ** 



** We abbreviate \upperJbound(cf)~\.^ by Ni. ** 



** The values in the matrix are calculated by multiplying the clock ** 
** probabilities by a factor of p, where p is the probability of ** 
** entering the state, and adding this value to the value already in ** 
** the position. ** 

procedure new state -matrix (matrix (cs, ct),p) 
begin procedure 
VI < h < Ni 



VI < i n < N n .matrix(cs,ct)[ii, . . . ,i n ] := 

n 

matrix(cs,ct)[ii, ...,i n ] + (px Y\p r ( c t e [** ~ <M*)) ) 

i=i 

end procedure 



